import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.spatial import distance
from tqdm import trange
import datetime


plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
plt.rc('font', size=10)
plt.rc('font', family='SimHei')

df = pd.read_excel("附件.xlsx")
x = df['x坐标 (m)'].values
y = df['y坐标 (m)'].values
vernal_equinox = datetime.date(2023, 3, 21)  # 2023年的春分

def find_nearest_neighbors(mirror_coords, k=6):
    """
    计算每个镜子最近的k个邻居

    参数:
        mirror_coords: 镜子坐标数组，形状为(n, 2)或(n, 3)
        k: 要查找的最近邻居数量

    返回:
        包含每个镜子最近k个邻居索引和距离的字典
    """
    # 计算所有镜子之间的欧式距离矩阵
    dist_matrix = distance.cdist(mirror_coords, mirror_coords, 'euclidean')

    # 将对角线设置为无穷大，避免每个镜子选择自己作为最近邻居
    np.fill_diagonal(dist_matrix, np.inf)

    # 获取每个镜子最近的k个邻居的索引和距离
    nearest_neighbors = []
    for i in range(len(mirror_coords)):
        # 获取距离最近的k个邻居的索引
        indices = np.argpartition(dist_matrix[i], k)[:k]
        # 获取对应的距离
        distances = dist_matrix[i][indices]
        # 按距离排序
        sorted_indices = indices[np.argsort(distances)]
        # 存储结果
        nearest_neighbors.append(sorted_indices.tolist())

    return nearest_neighbors

mirror_coords = df[['x坐标 (m)', 'y坐标 (m)']].values
neighbors = find_nearest_neighbors(mirror_coords, k=6)  # 每一个镜子他周围的最近六个镜子

phi = 39.4 / 360 * 2 * np.pi  # 北纬，单位为rad
H = 3000  # 海拔，单位为m
r = 350  # 半径，单位为m
h_tower = 80  # 塔的高度，单位为m
h_heater = 8  # 集热器高度，单位为m
D = 7  # 集热器直径，单位为m
r_no = 100  # 不安装定日镜的范围，单位为m
side_range = [2, 8]  # 定日镜边长在2m至8m之间
height_range = [2, 6]  # 安装高度在2m至6m之间
least_interval = 5  # 相邻定日镜底座中心之间的距离比镜面宽度多5m以上
G_0 = 1.366  # 太阳常数，单位为kW/m^2
N = len(df)  # 定日镜的总面数
n = 3  # 将每个定日镜反射的光认为是集中在它上面的n*n个点射出
sigma = 4.65e-3
tau = [0, 2 * np.pi / 5, 4 * np.pi / 5, 6 * np.pi / 5, 8 * np.pi / 5]
tau_len = len(tau)
time_list = [9, 10.5, 12, 13.5, 15]


# 太阳时角，ST为当地时间
def omiga(ST):
    return np.pi / 12 * (ST - 12)


# 太阳赤纬角。其中day为以春分作为第0天起算的天数，例如，若春分是3月21日，则4月1日对应day=11。
def delta(day):
    sin_delta = round(np.sin(2 * np.pi * day / 365) * np.sin(2 * np.pi * 23.45 / 360), 5)
    return np.arcsin(sin_delta)


# 太阳高度角。phi为当地纬度
def alpha_s(day, ST):
    global phi
    sin_alpha = round(np.cos(delta(day)) * np.cos(phi) * np.cos(omiga(ST)) + np.sin(delta(day)) * np.sin(phi), 5)
    return np.arcsin(sin_alpha)


# 太阳方位角
def gamma_s(day, ST):
    global phi
    cos_gamma = round((np.sin(delta(day)) - np.sin(alpha_s(day, ST)) * np.sin(phi)) / (np.cos(alpha_s(day, ST)) * np.cos(phi)), 5)
    return np.arccos(cos_gamma)


# 在直角坐标系下太阳的方向向量
def sun_direction(day, ST):
    return np.array([[np.cos(alpha_s(day, ST)) * np.sin(gamma_s(day, ST)), np.cos(alpha_s(day, ST)) * np.cos(gamma_s(day, ST)), np.sin(alpha_s(day, ST))]])


# 计算第i个定日镜的方位角和俯仰角。h_mirror为该面镜子的高度
def angle(day, ST, h_mirror, i):
    centor_direction = np.array([-df['x坐标 (m)'][i], -df['y坐标 (m)'][i], h_tower + h_heater / 2 - h_mirror])
    centor_direction = centor_direction / np.linalg.norm(centor_direction)  # 由镜面中心指向集热器中心的向量
    normal_vector = sun_direction(day, ST)[0] + centor_direction
    normal_vector = normal_vector / np.linalg.norm(normal_vector)  # 镜面法向量
    return np.array([np.arctan2(normal_vector[0], normal_vector[1]), np.arcsin(normal_vector[2])])
dir_angle, pitch_angle = angle(12, 6, 4, 0)
normal_vector = np.array([np.cos(pitch_angle) * np.sin(dir_angle), np.cos(pitch_angle) * np.cos(dir_angle), np.sin(pitch_angle)])


# 法向直接辐射辐照度，单位为kW/m^2
def DNI(day, ST):
    global G_0, H
    a = 0.4237 - 0.00821 * ((6 - H / 1000) ** 2)
    b = 0.5055 + 0.00595 * ((6.5 - H / 1000) ** 2)
    c = 0.2711 + 0.01858 * ((2.5 - H / 1000) ** 2)
    return G_0 * (a + b * np.exp(-c / np.sin(alpha_s(day, ST))))


# 余弦效率
def eta_cos(day, ST, h_mirror, i):
    centor_direction = np.array([-df['x坐标 (m)'][i], -df['y坐标 (m)'][i], h_tower + h_heater / 2 - h_mirror])
    centor_direction = centor_direction / np.linalg.norm(centor_direction)  # 由镜面中心指向集热器中心的向量
    centor_vector_norm = np.linalg.norm(centor_direction)
    normal_vector = sun_direction(day, ST)[0] + centor_direction
    normal_vector = normal_vector / np.linalg.norm(normal_vector)  # 镜面法向量
    normal_vector_norm = np.linalg.norm(normal_vector)
    opposite_norm = np.linalg.norm(centor_direction - normal_vector)  # 两向量夹角对边的范数
    return (normal_vector_norm ** 2 + centor_vector_norm ** 2 - opposite_norm ** 2) / (2 * normal_vector_norm * centor_vector_norm)


# 大气透射率
def eta_at(h_mirror, i):
    centor_direction = np.array([-df['x坐标 (m)'][i], -df['y坐标 (m)'][i], h_tower + h_heater / 2 - h_mirror])  # 由镜面中心指向集热器中心的向量
    d_HR = np.linalg.norm(centor_direction)  # 镜面中心到集热器中心的距离
    return 0.99321 - 0.0001176 * d_HR + 1.97e-8 * (d_HR ** 2)


# 判断直线与轴线在 z 轴上的圆柱是否有交点。
#
# 参数：
# - o: 直线上的一个点（三维 numpy 数组，形状为 (3,)）。
# - d: 直线的方向向量（三维 numpy 数组，形状为 (3,)）。
#
# 返回：
# - 确实有交点就返回true，否则返回false

def line_cylinder_intersection_z_axis(o, d):
    # 提取 x 和 y 分量
    o_x, o_y, o_z = o[0][0], o[1][0], o[2][0]
    d_x, d_y, d_z = d[0][0], d[0][1], d[0][2]
    # print(o_x, o_y, o_z, d_x, d_y, d_z)
    r = 3.5
    h = 8
    z0 = 80  # 底面中心的 z 坐标
    # 二次方程的系数
    A = d_x**2 + d_y**2
    B = 2 * (o_x * d_x + o_y * d_y)
    C = o_x**2 + o_y**2 - r**2
    # 解二次方程
    discriminant = B**2 - 4 * A * C
    if discriminant < 0:
        return 0
    sqrt_discriminant = np.sqrt(discriminant)
    t1 = (-B + sqrt_discriminant) / (2 * A)
    t2 = (-B - sqrt_discriminant) / (2 * A)
    # print("___________________________________________________")
    # print(((o[0][0] + t1 * d[0][0]) ** 2 + (o[1][0] + t1 * d[0][1]) ** 2) ** 0.5)
    # print(((o[0][0] + t2 * d[0][0]) ** 2 + (o[1][0] + t2 * d[0][1]) ** 2) ** 0.5)
    #
    # print(o[2][0]+t1*d[0][2])
    # print(o[2][0]+t2*d[0][2])
    # 计算交点
    for t in [t1, t2]:
        # 检查 z 坐标是否在 [z0, z0 + h] 范围内
        if z0 <= o[2][0]+t*d[0][2] <= z0 + h:
            return 1
    return 0


# 将在地面坐标系下的太阳向量转换到镜高为height、镜宽为width的第i面定日镜的镜面坐标上。E_H代表俯仰角
# 计算第i个定日镜的方位角和俯仰角。h_mirror为该面镜子的高度。angle(day, ST, h_mirror, i)
def sun_to_mirror(day, ST, h_mirror, i, sun_label):
    global n
    A_H, E_H = angle(day, ST, h_mirror, i)
    T = np.array([
        [-np.sin(E_H), -np.sin(A_H) * np.cos(E_H), np.cos(A_H) * np.cos(E_H)],
        [np.cos(E_H), -np.sin(A_H) * np.sin(E_H), np.cos(A_H) * np.sin(E_H)],
        [0, np.cos(A_H), np.sin(A_H)]
    ]).T
    return T @ sun_label.T


# 将镜高为height、镜宽为width的第i面定日镜上第j行、第k列的点的镜坐标转换为地面坐标。E_H代表俯仰角
# 计算第i个定日镜的方位角和俯仰角。h_mirror为该面镜子的高度。angle(day, ST, h_mirror, i)
def mirror_to_ground(day, ST, height, width, h_mirror, i, j, k):
    global n
    A_H, E_H = angle(day, ST, h_mirror, i)
    T = np.array([
        [-np.sin(E_H), -np.sin(A_H) * np.cos(E_H), np.cos(A_H) * np.cos(E_H)],
        [np.cos(E_H), -np.sin(A_H) * np.sin(E_H), np.cos(A_H) * np.sin(E_H)],
        [0, np.cos(A_H), np.sin(A_H)]
    ])
    mirror_label = np.array([[width / n * k + width / n * 1 / 2 - width / 2, height / n * j + height / n * 1 / 2 - height / 2, 0]]).T
    return T @ mirror_label + np.array([[df['x坐标 (m)'][i], df['y坐标 (m)'][i], h_mirror]]).T
# line_cylinder_intersection_z_axis(mirror_to_ground(31, 9, 6, 6, 4, 0, 0, 0), sun_direction(31, 9))
# count = 0
# for i in range(N):
#     for j in range(n):
#         for k in range(n):
#             count += line_cylinder_intersection_z_axis(mirror_to_ground(310, 9, 6, 6, 4, i, j, k), sun_direction(310, 9))
#             # count += line_cylinder_intersection_z_axis(np.array([[df['x坐标 (m)'][i], df['y坐标 (m)'][i], 4]]).T,
#             #                                            sun_direction(31, 9)[0], [0, 0, h_tower], D / 2, h_heater)
# print(count)


# 将在地面坐标系下坐标为groud_label的点转换到镜高为height、镜宽为width的第i面定日镜的镜面坐标上。E_H代表俯仰角
# 计算第i个定日镜的方位角和俯仰角。h_mirror为该面镜子的高度。angle(day, ST, h_mirror, i)
def ground_to_mirror(day, ST, h_mirror, i, ground_label):
    global n
    A_H, E_H = angle(day, ST, h_mirror, i)
    T = np.array([
        [-np.sin(E_H), -np.sin(A_H) * np.cos(E_H), np.cos(A_H) * np.cos(E_H)],
        [np.cos(E_H), -np.sin(A_H) * np.sin(E_H), np.cos(A_H) * np.sin(E_H)],
        [0, np.cos(A_H), np.sin(A_H)]
    ]).T
    return T @ (ground_label - np.array([[df['x坐标 (m)'][i], df['y坐标 (m)'][i], h_mirror]]).T)


# 判断阳光入射到第i面定日镜上第j行、第k列的点上时是否会被第m面镜子遮挡。
def occlusion_detection_incident(day, ST, height_i, width_i, height_m, width_m, h_mirror_i, h_mirror_m, i, j, k, m):
    ground_label = mirror_to_ground(day, ST, height_i, width_i, h_mirror_i, i, j, k)  # 该点在地面坐标系下的坐标
    mirror_label = ground_to_mirror(day, ST, h_mirror_m, m, ground_label)  # 该点在第m面定日镜的镜面坐标系下的坐标
    sun_direction_mirror = sun_to_mirror(day, ST, h_mirror_m, m, sun_direction(day, ST))  # 太阳光线的方向向量在第m面定日镜的镜面坐标系下的坐标
    x_2 = (sun_direction_mirror[2, 0] * mirror_label[0, 0] - sun_direction_mirror[0, 0] * mirror_label[2, 0]) / \
          sun_direction_mirror[2, 0]
    y_2 = (sun_direction_mirror[2, 0] * mirror_label[1, 0] - sun_direction_mirror[1, 0] * mirror_label[2, 0]) / \
          sun_direction_mirror[2, 0]
    if x_2 <= width_m / 2 and x_2 >= -width_m / 2 and y_2 <= height_m / 2 and y_2 >= -height_m / 2:
        return 1
    return 0


# 判断阳光从第i面定日镜上第j行、第k列的点出射时是否会被第m面镜子遮挡。
def occlusion_detection_emergent(day, ST, height_i, width_i, height_m, width_m, h_mirror_i, h_mirror_m, i, j, k, m):
    ground_label = mirror_to_ground(day, ST, height_i, width_i, h_mirror_i, i, j, k)  # 该点在地面坐标系下的坐标
    mirror_label = ground_to_mirror(day, ST, h_mirror_m, m, ground_label)  # 该点在第m面定日镜的镜面坐标系下的坐标

    centor_direction = np.array([-df['x坐标 (m)'][i], -df['y坐标 (m)'][i], h_tower + h_heater / 2 - h_mirror_i])
    centor_direction = np.array([centor_direction / np.linalg.norm(centor_direction)]).T  # 由镜面中心指向集热器中心的向量
    centor_direction = sun_to_mirror(day, ST, h_mirror_m, m, centor_direction.T)
    x_2 = (centor_direction[2, 0] * mirror_label[0, 0] - centor_direction[0, 0] * mirror_label[2, 0]) / \
          centor_direction[2, 0]
    y_2 = (centor_direction[2, 0] * mirror_label[1, 0] - centor_direction[1, 0] * mirror_label[2, 0]) / \
          centor_direction[2, 0]
    if x_2 <= width_m / 2 and x_2 >= -width_m / 2 and y_2 <= height_m / 2 and y_2 >= -height_m / 2:
        return 1
    return 0


def cone_to_ground(i, h_mirror_i, cone_label):
    centor_direction = np.array([-df['x坐标 (m)'][i], -df['y坐标 (m)'][i], h_tower + h_heater / 2 - h_mirror_i])
    centor_direction = np.array([centor_direction / np.linalg.norm(centor_direction)]).T# 由镜面中心指向集热器中心的向量

    ah, eh = np.array([np.arctan2(centor_direction[0], centor_direction[1]), np.arcsin(centor_direction[2])])
    ah = ah[0]
    eh = eh[0]
    # print(ah,eh)
    T = np.array([
            [-np.cos(ah), -np.sin(eh) * np.sin(ah), np.cos(eh) * np.sin(ah)],
            [np.sin(ah), -np.sin(eh) * np.cos(ah), np.cos(eh) * np.cos(ah)],
            [0, np.cos(eh), np.sin(eh)]
        ])
    return T @ cone_label


# 判断阳光从第i面定日镜上第j行、第k列的点出射后，发散为圆锥，以上面的四根光线为准，这些光线中的第a个是否可以照射到集热器上
def truncation_detection(day, ST, height_i, width_i, h_mirror_i, i, j, k, a):
    global sigma, tau

    ground_label = mirror_to_ground(day, ST, height_i, width_i, h_mirror_i, i, j, k)  # 该点在地面坐标系下的坐标

    cone_label = np.array([[np.sin(sigma) * np.cos(tau[a]), np.sin(sigma) * np.sin(tau[a]), np.cos(sigma)]]).T  # 圆锥坐标系下光线的向量
    # print('cone_label1', cone_label)
    cone_label = cone_to_ground(i, h_mirror_i, cone_label).T  # 将光线由圆锥坐标系转换到地面坐标系后的向量
    # print('cone_label2', cone_label)
    centor_direction = np.array([-df['x坐标 (m)'][i], -df['y坐标 (m)'][i], h_tower + h_heater / 2 - 4])

    centor_direction = np.array([centor_direction / np.linalg.norm(centor_direction)]).T  # 由镜面中心指向集热器中心的向量
    # print('centor_direction', centor_direction)

    return line_cylinder_intersection_z_axis(ground_label, cone_label)
# count = 0
# for i in range(N):
#     for j in range(n):
#         for k in range(n):
#             for a in range(4):
#                 count += truncation_detection(31, 12, 6, 6, 4, i, j, k, a)
# print(count)
# count = 0
# for i in range(N):
#     centor_direction = np.array([-df['x坐标 (m)'][i], -df['y坐标 (m)'][i], h_tower + h_heater / 2 - 4])
#     centor_direction = np.array([centor_direction / np.linalg.norm(centor_direction)]).T  # 由镜面中心指向集热器中心的向量
#     count += line_cylinder_intersection_z_axis(np.array([[df['x坐标 (m)'][i], df['y坐标 (m)'][i], 4]]).T, centor_direction.T)
# print(count)

# 第i面定日镜上第j行、第k列个点的光学效率
def eta(day, ST, h_mirror, i, j, k):
    eta_sb = 0  # 阴影遮挡效率
    flag = True
    ground_label = mirror_to_ground(day, ST, 6, 6, 4, i, j, k)  # 该点在地面坐标系下的坐标
    if line_cylinder_intersection_z_axis(ground_label, sun_direction(day, ST)) != 1:
        for m in range(6):
            if occlusion_detection_incident(day, ST, 6, 6, 6, 6, 4, 4, i, j, k, neighbors[i][m]) == 1:
                flag = False
                break
            elif occlusion_detection_emergent(day, ST, 6, 6, 6, 6, 4, 4, i, j, k, neighbors[i][m]) == 1:
                flag = False
    if flag:
        eta_sb = 1

    count = 0
    for a in range(tau_len):
        count += truncation_detection(day, ST, 6, 6, 4, i, j, k, a)

    eta_trunc = count / tau_len  # 集热器截断效率
    eta_ref = 0.92  # 镜面反射率，可取为常数
    return eta_sb * eta_cos(day, ST, h_mirror, i) * eta_at(h_mirror, i) * eta_trunc * eta_ref, eta_cos(day, ST, h_mirror, i), eta_sb, eta_trunc


# 定日镜采光面积。单位为m^2
def A():
    return (6 / n) * (6 / n)


# 定日镜场的输出热功率
def E_field():
    global N, n
    eta_total_value = []
    eta_cos_total_value = []
    eta_sb_total_value = []
    eta_trunc_total_value = []
    for i in trange(N):
        eta_total = 0
        eta_cos_total = 0
        eta_sb_total = 0
        eta_trunc_total = 0
        for m in range(12):
            day = (datetime.date(2023, m + 1, 21) - vernal_equinox).days
            for ST in time_list:
                for j in range(n):
                    for k in range(n):
                        delta_eta = eta(day, ST, 4, i, j, k)
                        eta_total += delta_eta[0]
                        eta_cos_total += delta_eta[1]
                        eta_sb_total += delta_eta[2]
                        eta_trunc_total += delta_eta[3]

        eta_total_value.append(eta_total / n / n / 5 / 12)
        eta_cos_total_value.append(eta_cos_total / n / n / 5 / 12)
        eta_sb_total_value.append(eta_sb_total / n / n / 5 / 12)
        eta_trunc_total_value.append(eta_trunc_total / n / n / 5 / 12)

    # 创建图形
    plt.figure(figsize=(8, 6))
    # 绘制散点图，使用颜色映射
    scatter = plt.scatter(x, y, c=eta_total_value, cmap='viridis', s=100, alpha=0.8)
    # 添加颜色条
    plt.colorbar(scatter, label='值')
    # 添加标题和轴标签
    plt.title('总光学效率分布图')
    plt.xlabel('X坐标')
    plt.ylabel('Y坐标')
    # 显示图形
    plt.show()

    # 创建图形
    plt.figure(figsize=(8, 6))
    # 绘制散点图，使用颜色映射
    scatter = plt.scatter(x, y, c=eta_cos_total_value, cmap='viridis', s=100, alpha=0.8)
    # 添加颜色条
    plt.colorbar(scatter, label='值')
    # 添加标题和轴标签
    plt.title('总余弦效率分布图')
    plt.xlabel('X坐标')
    plt.ylabel('Y坐标')
    # 显示图形
    plt.show()

    # 创建图形
    plt.figure(figsize=(8, 6))
    # 绘制散点图，使用颜色映射
    scatter = plt.scatter(x, y, c=eta_sb_total_value, cmap='viridis', s=100, alpha=0.8)
    # 添加颜色条
    plt.colorbar(scatter, label='值')
    # 添加标题和轴标签
    plt.title('阴影遮挡效率分布图')
    plt.xlabel('X坐标')
    plt.ylabel('Y坐标')
    # 显示图形
    plt.show()

    # 创建图形
    plt.figure(figsize=(8, 6))
    # 绘制散点图，使用颜色映射
    scatter = plt.scatter(x, y, c=eta_trunc_total_value, cmap='viridis', s=100, alpha=0.8)
    # 添加颜色条
    plt.colorbar(scatter, label='值')
    # 添加标题和轴标签
    plt.title('集热器截断效率分布图')
    plt.xlabel('X坐标')
    plt.ylabel('Y坐标')
    # 显示图形
    plt.show()

E_field()
